A high-throughput phenotyping pipeline for quinoa (Chenopodium quinoa) panicles using image analysis with convolutional neural networks

1 Setup

source('https://inkaverse.com/setup.r')
library(readxl)

session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.0 (2024-04-24 ucrt)
 os       Windows 11 x64 (build 22631)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  Spanish_Latin America.utf8
 ctype    Spanish_Latin America.utf8
 tz       America/Lima
 date     2024-06-16
 pandoc   3.1.11 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package       * version  date (UTC) lib source
 agricolae       1.3-7    2023-10-22 [1] CRAN (R 4.4.0)
 AlgDesign       1.2.1    2022-05-25 [1] CRAN (R 4.4.0)
 askpass         1.2.0    2023-09-03 [1] CRAN (R 4.4.0)
 boot            1.3-30   2024-02-26 [2] CRAN (R 4.4.0)
 cachem          1.1.0    2024-05-16 [1] CRAN (R 4.4.0)
 cellranger      1.1.0    2016-07-27 [1] CRAN (R 4.4.0)
 cli             3.6.2    2023-12-11 [1] CRAN (R 4.4.0)
 cluster         2.1.6    2023-12-01 [2] CRAN (R 4.4.0)
 codetools       0.2-20   2024-03-31 [2] CRAN (R 4.4.0)
 colorspace      2.1-0    2023-01-23 [1] CRAN (R 4.4.0)
 cowplot       * 1.1.3    2024-01-22 [1] CRAN (R 4.4.0)
 curl            5.2.1    2024-03-01 [1] CRAN (R 4.4.0)
 devtools      * 2.4.5    2022-10-11 [1] CRAN (R 4.4.0)
 digest          0.6.35   2024-03-11 [1] CRAN (R 4.4.0)
 dplyr         * 1.1.4    2023-11-17 [1] CRAN (R 4.4.0)
 DT              0.33     2024-04-04 [1] CRAN (R 4.4.0)
 ellipsis        0.3.2    2021-04-29 [1] CRAN (R 4.4.0)
 emmeans         1.10.2   2024-05-20 [1] CRAN (R 4.4.0)
 estimability    1.5.1    2024-05-12 [1] CRAN (R 4.4.0)
 evaluate        0.24.0   2024-06-10 [1] CRAN (R 4.4.0)
 FactoMineR    * 2.11     2024-04-20 [1] CRAN (R 4.4.0)
 fansi           1.0.6    2023-12-08 [1] CRAN (R 4.4.0)
 fastmap         1.2.0    2024-05-15 [1] CRAN (R 4.4.0)
 flashClust      1.01-2   2012-08-21 [1] CRAN (R 4.4.0)
 forcats       * 1.0.0    2023-01-29 [1] CRAN (R 4.4.0)
 fs              1.6.4    2024-04-25 [1] CRAN (R 4.4.0)
 gargle          1.5.2    2023-07-20 [1] CRAN (R 4.4.0)
 generics        0.1.3    2022-07-05 [1] CRAN (R 4.4.0)
 ggplot2       * 3.5.1    2024-04-23 [1] CRAN (R 4.4.0)
 ggrepel         0.9.5    2024-01-10 [1] CRAN (R 4.4.0)
 glue            1.7.0    2024-01-09 [1] CRAN (R 4.4.0)
 googledrive   * 2.1.1    2023-06-11 [1] CRAN (R 4.4.0)
 googlesheets4 * 1.1.1    2023-06-11 [1] CRAN (R 4.4.0)
 gsheet        * 0.4.5    2020-04-07 [1] CRAN (R 4.4.0)
 gtable          0.3.5    2024-04-22 [1] CRAN (R 4.4.0)
 hms             1.1.3    2023-03-21 [1] CRAN (R 4.4.0)
 htmltools       0.5.8.1  2024-04-04 [1] CRAN (R 4.4.0)
 htmlwidgets     1.6.4    2023-12-06 [1] CRAN (R 4.4.0)
 httpuv          1.6.15   2024-03-26 [1] CRAN (R 4.4.0)
 httr            1.4.7    2023-08-15 [1] CRAN (R 4.4.0)
 huito         * 0.2.4    2023-10-25 [1] CRAN (R 4.4.0)
 inti          * 0.6.5    2024-05-16 [1] CRAN (R 4.4.0)
 jsonlite        1.8.8    2023-12-04 [1] CRAN (R 4.4.0)
 knitr         * 1.47     2024-05-29 [1] CRAN (R 4.4.0)
 later           1.3.2    2023-12-06 [1] CRAN (R 4.4.0)
 lattice         0.22-6   2024-03-20 [2] CRAN (R 4.4.0)
 leaps           3.2      2024-06-10 [1] CRAN (R 4.4.0)
 lifecycle       1.0.4    2023-11-07 [1] CRAN (R 4.4.0)
 lme4            1.1-35.3 2024-04-16 [1] CRAN (R 4.4.0)
 lubridate     * 1.9.3    2023-09-27 [1] CRAN (R 4.4.0)
 magick        * 2.8.3    2024-02-18 [1] CRAN (R 4.4.0)
 magrittr        2.0.3    2022-03-30 [1] CRAN (R 4.4.0)
 MASS            7.3-60.2 2024-04-24 [2] local
 Matrix          1.7-0    2024-03-22 [2] CRAN (R 4.4.0)
 memoise         2.0.1    2021-11-26 [1] CRAN (R 4.4.0)
 mime            0.12     2021-09-28 [1] CRAN (R 4.4.0)
 miniUI          0.1.1.1  2018-05-18 [1] CRAN (R 4.4.0)
 minqa           1.2.7    2024-05-20 [1] CRAN (R 4.4.0)
 mnormt          2.1.1    2022-09-26 [1] CRAN (R 4.4.0)
 multcomp        1.4-25   2023-06-20 [1] CRAN (R 4.4.0)
 multcompView    0.1-10   2024-03-08 [1] CRAN (R 4.4.0)
 munsell         0.5.1    2024-04-01 [1] CRAN (R 4.4.0)
 mvtnorm         1.2-5    2024-05-21 [1] CRAN (R 4.4.0)
 nlme            3.1-164  2023-11-27 [2] CRAN (R 4.4.0)
 nloptr          2.0.3    2022-05-26 [1] CRAN (R 4.4.0)
 openssl         2.2.0    2024-05-16 [1] CRAN (R 4.4.0)
 pillar          1.9.0    2023-03-22 [1] CRAN (R 4.4.0)
 pkgbuild        1.4.4    2024-03-17 [1] CRAN (R 4.4.0)
 pkgconfig       2.0.3    2019-09-22 [1] CRAN (R 4.4.0)
 pkgload         1.3.4    2024-01-16 [1] CRAN (R 4.4.0)
 profvis         0.3.8    2023-05-02 [1] CRAN (R 4.4.0)
 promises        1.3.0    2024-04-05 [1] CRAN (R 4.4.0)
 psych         * 2.4.3    2024-03-18 [1] CRAN (R 4.4.0)
 purrr         * 1.0.2    2023-08-10 [1] CRAN (R 4.4.0)
 R6              2.5.1    2021-08-19 [1] CRAN (R 4.4.0)
 rappdirs        0.3.3    2021-01-31 [1] CRAN (R 4.4.0)
 Rcpp            1.0.12   2024-01-09 [1] CRAN (R 4.4.0)
 readr         * 2.1.5    2024-01-10 [1] CRAN (R 4.4.0)
 readxl        * 1.4.3    2023-07-06 [1] CRAN (R 4.4.0)
 remotes         2.5.0    2024-03-17 [1] CRAN (R 4.4.0)
 rlang           1.1.4    2024-06-04 [1] CRAN (R 4.4.0)
 rmarkdown       2.27     2024-05-17 [1] CRAN (R 4.4.0)
 rstudioapi      0.16.0   2024-03-24 [1] CRAN (R 4.4.0)
 sandwich        3.1-0    2023-12-11 [1] CRAN (R 4.4.0)
 scales          1.3.0    2023-11-28 [1] CRAN (R 4.4.0)
 scatterplot3d   0.3-44   2023-05-05 [1] CRAN (R 4.4.0)
 sessioninfo     1.2.2    2021-12-06 [1] CRAN (R 4.4.0)
 shiny         * 1.8.1.1  2024-04-02 [1] CRAN (R 4.4.0)
 showtext        0.9-7    2024-03-02 [1] CRAN (R 4.4.0)
 showtextdb      3.0      2020-06-04 [1] CRAN (R 4.4.0)
 stringi         1.8.4    2024-05-06 [1] CRAN (R 4.4.0)
 stringr       * 1.5.1    2023-11-14 [1] CRAN (R 4.4.0)
 survival        3.5-8    2024-02-14 [2] CRAN (R 4.4.0)
 sysfonts        0.8.9    2024-03-02 [1] CRAN (R 4.4.0)
 TH.data         1.1-2    2023-04-17 [1] CRAN (R 4.4.0)
 tibble        * 3.2.1    2023-03-20 [1] CRAN (R 4.4.0)
 tidyr         * 1.3.1    2024-01-24 [1] CRAN (R 4.4.0)
 tidyselect      1.2.1    2024-03-11 [1] CRAN (R 4.4.0)
 tidyverse     * 2.0.0    2023-02-22 [1] CRAN (R 4.4.0)
 timechange      0.3.0    2024-01-18 [1] CRAN (R 4.4.0)
 tzdb            0.4.0    2023-05-12 [1] CRAN (R 4.4.0)
 urlchecker      1.0.1    2021-11-30 [1] CRAN (R 4.4.0)
 usethis       * 2.2.3    2024-02-19 [1] CRAN (R 4.4.0)
 utf8            1.2.4    2023-10-22 [1] CRAN (R 4.4.0)
 vctrs           0.6.5    2023-12-01 [1] CRAN (R 4.4.0)
 withr           3.0.0    2024-01-16 [1] CRAN (R 4.4.0)
 xfun            0.44     2024-05-15 [1] CRAN (R 4.4.0)
 xtable          1.8-4    2019-04-21 [1] CRAN (R 4.4.0)
 yaml            2.3.8    2023-12-11 [1] CRAN (R 4.4.0)
 zoo             1.8-12   2023-04-13 [1] CRAN (R 4.4.0)

 [1] C:/Users/floza/AppData/Local/R/win-library/4.4
 [2] C:/Program Files/R/R-4.4.0/library

──────────────────────────────────────────────────────────────────────────────

2 Googlesheets connect

url <- "https://docs.google.com/spreadsheets/d/1ywk2jhCKRKy_9-8pH9LPjBN8jmfMr_6ofBB7NY0WY5Y"
# browseURL(gsh)
gs <- as_sheets_id(url)

3 Import pipeline results

Import result from the serve

3.1 Classification

ssh flaaw57@server231.ipsp.uni-hohenheim.de
cd /work/workspaces/flavio/segmentation/img-segmented/1-segmentation-part

cp -r /work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/ /work/workspaces/flavio/segmentation/img-segmented

Import result to local

cd /work/workspaces/flavio/segmentation/img-segmented/1-segmentation-part/

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/flavio/segmentation/img-segmented/1-segmentation-part/classification_prediction.csv 2-stage-approach/2-classification-part/

rs_cla <- list.files(path = "classification", pattern = "classification_prediction.csv"
                     , recursive = T, full.names = T) %>% 
  read_csv() %>% 
  mutate(info = gsub("^.*part\\/", "", Imagen)) %>% 
  separate(info, c("img_folder", "type", "image"), sep = "/") %>% 
  separate(image, c("img_name", "panicle_number"), sep = "_panicle_") %>% 
  mutate(panicle_number = gsub("\\D", "", panicle_number)) %>% 
  separate(img_folder, c("img_folder", "nfolder"), sep = "_")

rs_cla$img_folder %>% unique() %>% sort()
## [1] "F6-camacani"  "F7-illpa"     "F8-camacani"  "F8-illpa"     "scale-images"

3.1.1 Classiciation pipeline

• Server: quino-img-analysis/2-stage-approach/2-classification-part/classification_quinoa-panicles.ipynb
• Local: classification/quinoa-panicles_classification-model_accuracy.ipynb

3.2 Segementation

Import result from the serve

ssh flaaw57@server231.ipsp.uni-hohenheim.de

find /work/workspaces/flavio/segmentation/results/ -name "*.csv" -exec cp {} /work/workspaces/flavio/segmentation/results/csvresults/ \;

cd /work/workspaces/flavio/segmentation/results/csvresults

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/flavio/segmentation/results/csvresults/ 2-stage-approach/1-segmentation-part/result-final/

ssh flaaw57@server231.ipsp.uni-hohenheim.de
cd /work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/F7-illpa_3 xclude

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/F7-illpa_4 xclude

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/F7-illpa_5 xclude

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/F8-camacani_1 xclude

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/F8-illpa xclude

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/quinoa-image-analysis/2-quinoa-panicle-maskrcnn/0-quinoa-panicle-stacked/1-segmentation-part/scale-images_1 xclude

• Copy results in csv to git project

Manual rename folder scale-images_1 > F9-camamacani_1

seg.files <- list.files("xclude/", pattern = ".csv"
                        , recursive = T, full.names = T)

cp <- 1:length(seg.files) %>% 
  map(\(x) {
    
    opath <- seg.files[x]
    
    dir <-  gsub(".*/([^/]+)/.*", "\\1", opath)
    
    fname <- basename(opath)
    
    file.copy(from = opath
              , to = file.path("2-stage-approach/1-segmentation-part/results/"
                               , paste0(dir, "-", fname)))
  })

seg.files <- list.files("segmentation/results/", pattern = ".*result.*.csv"
                        , recursive = T, full.names = T)

rs_seg <- 1:length(seg.files) %>% map( \(x) {
  
  seg.files[x] %>% 
    read.csv()
  
}) %>% 
  bind_rows() %>% 
  separate(img_folder, c("img_folder", "nfolder"), sep = "_")
  
rs_seg$img_folder %>% unique() %>% sort()
## [1] "F6-camacani"  "F7-illpa"     "F8-camacani"  "F8-illpa"     "scale-images"

3.2.1 Classiciation pipeline

• Lydia (missing)

3.3 Merge data: segmentation + classification

rs_cla %>% names()
## [1] "Imagen"         "Accuracy"       "Clase"          "img_folder"    
## [5] "nfolder"        "type"           "img_name"       "panicle_number"
rs_seg %>% names()
##  [1] "img_folder"        "nfolder"           "img_name"         
##  [4] "img_qrcode"        "img_res"           "scale_pixelpercm2"
##  [7] "scale_RGB_R"       "scale_RGB_G"       "scale_RGB_B"      
## [10] "scale_amount"      "scale_median"      "scale_mean"       
## [13] "scale_std"         "scale_min"         "scale_max"        
## [16] "panicle_detected"  "panicle_number"    "panicle_width"    
## [19] "panicle_length"    "panicle_area"      "panicle_shape"    
## [22] "panicle_RGB_mean"  "panicle_RGB_stdev"

rs <- rs_cla %>% 
  merge(rs_seg, by = c("img_folder", "img_name", "panicle_number")) %>% 
  dplyr::select(!c(panicle_shape, starts_with("nfolder"))) %>% 
  rename(panicle_shape = "Clase") %>% 
  mutate(img_folder = gsub("\\_[0-9]+$", "", img_folder)) %>% 
  separate(panicle_RGB_mean, c("panicle_R_mean", "panicle_G_mean", "panicle_B_mean"), sep = ":") %>% 
  separate(panicle_RGB_stdev, c("panicle_R_stdev", "panicle_G_stdev", "panicle_B_stdev"), sep = ":") %>% 
  mutate(`panicle_length.width` = panicle_length/panicle_width) %>% 
  mutate(`panicle_width.length` = panicle_width/panicle_length)

rs %>% names()
##  [1] "img_folder"           "img_name"             "panicle_number"      
##  [4] "Imagen"               "Accuracy"             "panicle_shape"       
##  [7] "type"                 "img_qrcode"           "img_res"             
## [10] "scale_pixelpercm2"    "scale_RGB_R"          "scale_RGB_G"         
## [13] "scale_RGB_B"          "scale_amount"         "scale_median"        
## [16] "scale_mean"           "scale_std"            "scale_min"           
## [19] "scale_max"            "panicle_detected"     "panicle_width"       
## [22] "panicle_length"       "panicle_area"         "panicle_R_mean"      
## [25] "panicle_G_mean"       "panicle_B_mean"       "panicle_R_stdev"     
## [28] "panicle_G_stdev"      "panicle_B_stdev"      "panicle_length.width"
## [31] "panicle_width.length"
rs$img_folder %>% unique()
## [1] "F6-camacani"  "F7-illpa"     "F8-camacani"  "F8-illpa"     "scale-images"

rs %>% 
  web_table(caption = "Classification and Segementation results")

3.4 Model classification accurary

ssh flaaw57@server231.ipsp.uni-hohenheim.de
cd /work/workspaces/flavio/classification
cd /work/workspaces/flavio/paniclesimgs

Activate Remote jupyeter notebook

screen
jupyter notebook --no-browser --port=8080
# Detaching a Screen > ctrl + a + d

ssh -L 8080:localhost:8080 flaaw57@server231.ipsp.uni-hohenheim.de

Open: http://localhost:8080/

Check/Stop running process

screen -ls
screen -r 47624
# 
#> Detaching a Screen > ctrl + a + d

ps aux | grep jupyter
# kill 27247

Save the result as .txt file and download to git repository

ssh flaaw57@server231.ipsp.uni-hohenheim.de
cd /work/workspaces/flavio/classification

ls > models_classification.txt

scp -r flaaw57@server231.ipsp.uni-hohenheim.de:/work/workspaces/flavio/classification/models_classification.txt 2-stage-approach/2-classification-part/

mdl_class_info <- gs %>% 
  range_read("ModelClass") %>% 
  rename_with(~ tolower(gsub(".", "_", .x, fixed = TRUE)))

mdl_class <- list.files("classification"
                        , pattern = ".*models_classification.*.txt"
                        , recursive = T, full.names = T) %>% 
  read.delim(header = F) %>% 
  rename(files = V1) %>% 
  filter(str_detect(files, pattern = "best-model")) %>%
  mutate(across(everything(), ~gsub(".h5", "", .))) %>% 
  separate(files, c("architecture", "num", "epoch", "type", "accuracy"), sep = "_") %>% 
  merge(., mdl_class_info, by = c("architecture", "num"))

4 Manuscript

5 Table 1

sheet <- gs %>% 
  range_read("albums") 

tab <- sheet %>%
  filter(stage %in% "flowering" & type %in% "panicles") %>% 
  dplyr::select(contains(c("seas", "geno", "gene", "loc", "design", "dev", "res"))) %>% 
  add_column(Pictures = c(3862, 1240, 25, 24))

tab %>%
  web_table()

# tab %>% write_sheet(ss = gs, sheet = "tab1")

6 Figure 1

variation <- list.files(path = "annotation/variation"
                    , pattern = "jpg"
                    , recursive = T
                    , full.names = T
                    ) %>% 
  enframe(name = "img", value = "path") %>% 
  mutate(name = unglue::unglue_vec(path, "{}variation/{x}.jpg{}")) %>% 
  separate(name, c("name", "num"), sep = "_") %>% 
  mutate(group = case_when(
    str_detect(name, "damage|blurred|dry") ~ "excluded"
    , TRUE ~ "variation"
  )) %>% 
  mutate(plot = paste0("rasterGrob(image_read(","'", path,"'",") %>%  image_resize(., '400x550!'))")) %>% 
  arrange(desc(group)) %>% 
  mutate(name = gsub("-", " ", name))

rcd <- list("background brigtness" = "background brightness"
            , "presence of leaves" = "leaves within panicles"
            , "over dry panicles" = "overly dry images"
            )

var <- variation %>% 
  filter(group %in% "variation") %>% 
  mutate(across(name, ~ dplyr::recode(., !!!rcd))) %>% 
  mutate(across(name, ~ str_to_title(.)))
  

plot1 <- 1:length(unique(var$name)) %>% map(function(x) {
  
  names <- unique(var$name)[x]
  
  grid <- var %>% 
    filter(name %in% names) %>% 
    select(plot) %>% 
    deframe()
  
  list <- 1:length(grid) %>% map(function(y) {
    
     eval(parse(text = grid[y]))
    
  }) %>% 
    plot_grid(plotlist = .
            , ncol = 2
            , labels = c("", names)
            , label_x = -0.35
            , label_y = 1.02
            , label_size = 25
            , hjust = 0
            , vjust = 0
            ) +
    theme(plot.margin = unit(c(1, 0.2, 0, 0), "cm")) 
  
}) %>% 
  plot_grid(plotlist = . 
            , nrow =  2
            ) +
  theme(plot.margin = unit(c(1, 0, 0, 0), "cm")) 

var <- variation %>% 
  filter(group %in% "excluded") %>% 
  mutate(across(name, ~ dplyr::recode(., !!!rcd))) %>% 
  mutate(across(name, ~str_to_title(.)))

plot2 <- 1:length(unique(var$name)) %>% map(function(x) {
  
  names <- unique(var$name)[x]
  
  grid <- var %>% 
    filter(name %in% names) %>% 
    select(plot) %>% 
    deframe()
  
  list <- 1:length(grid) %>% map(function(y) {
    
     eval(parse(text = grid[y]))
    
  }) %>% 
    plot_grid(plotlist = .
            , ncol = 2
            , labels = c("", names)
            , label_x = -0.35
            , label_y = 1.02
            , label_size = 25
            , hjust = 0
            , vjust = 0
            ) +
     theme(plot.margin = unit(c(1, 0.2, 0, 0), "cm")) 
  
}) %>% 
  plot_grid(plotlist = . , nrow =  1) +
  theme(plot.margin = unit(c(1, 0, 0, 0), "cm")) 

plot <- plot_grid(plot1, plot2
          , nrow = 2
          , labels = c("(a) Variation among images", "(b) Excluded images")
          , rel_heights = c(2, 1)
          , hjust = -0.02
          , vjust = 1.1
          , label_size = 25
          ) 

plot %>% 
    ggsave2(filename = "files/Fig1.jpg"
          , plot = .
          , dpi = 200
          , width = 19.5
          , height = 15
          )

plot %>% 
    ggsave2(filename = "files/Fig1.eps"
          , plot = .
          , dpi = 200
          , width = 19.5
          , height = 15
          )

include_graphics("files/Fig1.jpg")

7 Table 2

tab <- "model-predictions/model-result-detections/models.xlsx" %>% 
  readxl::read_xlsx() %>%
  mutate(`mask-resolution` = case_when(
    `mask-resolution` == "64x64" ~ "56x56"
    , `mask-resolution` == "32x32" ~ "28x28"
    , TRUE ~ `mask-resolution`
  )) %>% 
  rename(Model = model
         , `imgaug` = `image-augmentation`
         , `mask resolution` = `mask-resolution`
         , `loss weight` = `loss-weight`
         , `heads.m` = layers
         ) 

tab %>%
  web_table()

# tab %>% write_sheet(ss = gs, sheet = "tab2")

8 Figure 2

8.1 ModelSeg

library(car)
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))

files <- list.files("segmentation/statistical-analysis/"
                    , pattern = "param.xlsx"
                    , full.names = T
                    , recursive = T
                    )

ap595 <- files %>% 
  readxl::read_xlsx(sheet = 2) %>% 
  mutate(across(1:binac_code, ~as.factor(.))) %>% 
  select(model = modeltype, loss, heads.m, mask_resolution, ap595) %>% 
  mutate(model = paste0("segmentation-", model))

# ap595 %>% write_sheet(gs, data = ., sheet = "stab1.0")

model <- ap595 %>% 
  lm(ap595 ~ loss*heads.m*mask_resolution, data = .)

Anova(model, type="III",singular.ok = TRUE)
## Anova Table (Type III tests)
## 
## Response: ap595
##                              Sum Sq Df    F value                Pr(>F)    
## (Intercept)                  44.991  1 60193.5877 < 0.00000000000000022 ***
## loss                          0.001  3     0.6188               0.60536    
## heads.m                       0.039  1    51.9138       0.0000000008114 ***
## mask_resolution               0.020  1    26.4429       0.0000027749151 ***
## loss:heads.m                  0.005  3     2.2900               0.08672 .  
## loss:mask_resolution          0.003  3     1.2449               0.30086    
## heads.m:mask_resolution       0.003  1     3.8630               0.05370 .  
## loss:heads.m:mask_resolution  0.009  3     4.0186               0.01104 *  
## Residuals                     0.048 64                                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
av <- anova(model)

tb <- Anova(model, type="III",singular.ok = TRUE) %>% 
  as.data.frame() %>% 
  rownames_to_column("Factor") %>% 
      mutate(Sig = case_when(
        `Pr(>F)` <= 0.001  ~ "***"
        , `Pr(>F)` <= 0.01  ~ "**"
        , `Pr(>F)` <= 0.05  ~ "*"
        , `Pr(>F)` > 0.05 ~ "ns"
      )) %>% 
  tibble() %>% 
  tibble::add_row(Factor = "---") %>%
  tibble::add_row(Factor = "Significance:") %>% 
  tibble::add_row(Factor = "0.001 '***' 0.01 '**' 0.05 '*'") 
  

# tb %>% write_sheet(ss = gs, sheet = "ModelSegAOV")

mc <- emmeans::emmeans(model, ~ loss | mask_resolution | heads.m ) %>% 
  multcomp::cld(Letters = letters, reversed = T) %>% 
  tibble()

# sqrt(tail(av$`Mean Sq`, 1))/mean(mc$emmean)*100

# mc %>% write_sheet(ss = gs, sheet = "ModelSeg")

plot1 <- mc %>% 
  plot_smr(type = "bar"
           , x = "loss"
           , xlab = "Loss weight" 
           , y = "emmean"
           , ylimits = c(0, 1.1, 0.1)
           , ylab = "mAP"
           , group = "mask_resolution"
           , glab = "Mask resolution"
           , sig = ".group"
           # , error = "SE"
           ) +
  facet_wrap(. ~ heads.m, ncol = 3) +
  theme(legend.position.inside = c(0.15, 0.97)
        , legend.direction = "horizontal"
        , strip.text = element_text(size = 12)
        , axis.text.x = element_text(size = 10)
        )

tab <- mc %>% 
  tibble() %>% 
  rename("loss weight" = "loss", "mask resolution" = "mask_resolution") %>% 
  mutate(Model = paste0("segmentation-", row_number()), .before = `loss weight`)

# tab %>% 
#   mutate(across(where(is.numeric), ~round(., 3))) %>% 
#   select(!c(df, lower.CL, upper.CL)) %>% 
#   rename(ste = "SE", sig = ".group", ap595 = "emmean") %>% 
#   mutate(across(sig, ~trimws(.))) %>% 
#   write_sheet(ss = gs, sheet = "tab2") 

8.2 ModelClass

library(car)
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))

model <- mdl_class %>% 
  lm(accuracy ~ architecture*`dense layers`*`activation function`, data = .)

# mdl_class %>% write_sheet(gs, data = ., sheet = "stab2.0")

Anova(model, type="III", singular.ok = TRUE)
## Anova Table (Type III tests)
## 
## Response: accuracy
##                                                    Sum Sq Df    F value
## (Intercept)                                       11.0405  1 15323.7045
## architecture                                       0.0012  2     0.8439
## `dense layers`                                     0.0001  1     0.0998
## `activation function`                              0.0012  1     1.6945
## architecture:`dense layers`                        0.0037  2     2.5392
## architecture:`activation function`                 0.0024  2     1.6638
## `dense layers`:`activation function`               0.0002  1     0.2472
## architecture:`dense layers`:`activation function`  0.0037  2     2.5410
## Residuals                                          0.0166 23           
##                                                                Pr(>F)    
## (Intercept)                                       <0.0000000000000002 ***
## architecture                                                   0.4429    
## `dense layers`                                                 0.7549    
## `activation function`                                          0.2059    
## architecture:`dense layers`                                    0.1008    
## architecture:`activation function`                             0.2114    
## `dense layers`:`activation function`                           0.6237    
## architecture:`dense layers`:`activation function`              0.1007    
## Residuals                                                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model)
## Analysis of Variance Table
## 
## Response: accuracy
##                                                   Df    Sum Sq    Mean Sq
## architecture                                       2 0.0019413 0.00097064
## `dense layers`                                     1 0.0000215 0.00002149
## `activation function`                              1 0.0023281 0.00232814
## architecture:`dense layers`                        2 0.0040591 0.00202953
## architecture:`activation function`                 2 0.0013703 0.00068516
## `dense layers`:`activation function`               1 0.0005054 0.00050538
## architecture:`dense layers`:`activation function`  2 0.0036616 0.00183078
## Residuals                                         23 0.0165712 0.00072049
##                                                   F value  Pr(>F)  
## architecture                                       1.3472 0.27972  
## `dense layers`                                     0.0298 0.86438  
## `activation function`                              3.2313 0.08538 .
## architecture:`dense layers`                        2.8169 0.08049 .
## architecture:`activation function`                 0.9510 0.40104  
## `dense layers`:`activation function`               0.7014 0.41092  
## architecture:`dense layers`:`activation function`  2.5410 0.10068  
## Residuals                                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tb <- Anova(model, type="III",singular.ok = TRUE) %>% 
  as.data.frame() %>% 
  rownames_to_column("Factor") %>% 
      mutate(Sig = case_when(
        `Pr(>F)` <= 0.001  ~ "***"
        , `Pr(>F)` <= 0.01  ~ "**"
        , `Pr(>F)` <= 0.05  ~ "*"
        , `Pr(>F)` > 0.05 ~ "ns"
      )) %>% 
  tibble() %>% 
  tibble::add_row(Factor = "---") %>%
  tibble::add_row(Factor = "Significance:") %>% 
  tibble::add_row(Factor = "0.001 '***' 0.01 '**' 0.05 '*'") 
  
# tb %>% write_sheet(ss = gs, sheet = "ModelClasAOV")

mc <- emmeans::emmeans(model, ~ architecture*`dense layers`*`activation function`) %>% 
  multcomp::cld(Letters = letters, reversed = T)

mc
##  architecture   dense layers activation function emmean     SE df lower.CL
##  InceptionV3            1024 softmax              0.924 0.0155 23    0.892
##  InceptionV3             128 sigmoid              0.917 0.0134 23    0.889
##  VGG16                  1024 softmax              0.909 0.0155 23    0.877
##  VGG16                   128 softmax              0.909 0.0155 23    0.877
##  InceptionV3            1024 sigmoid              0.909 0.0155 23    0.877
##  VGG16                   128 sigmoid              0.909 0.0155 23    0.877
##  EfficientNetB0          128 sigmoid              0.909 0.0155 23    0.877
##  EfficientNetB0          128 softmax              0.901 0.0190 23    0.862
##  VGG16                  1024 sigmoid              0.901 0.0190 23    0.862
##  EfficientNetB0         1024 sigmoid              0.901 0.0190 23    0.862
##  InceptionV3             128 softmax              0.871 0.0134 23    0.843
##  EfficientNetB0         1024 softmax              0.854 0.0155 23    0.821
##  upper.CL .group
##     0.956  a    
##     0.944  a    
##     0.941  a    
##     0.941  a    
##     0.941  a    
##     0.941  a    
##     0.941  a    
##     0.941  a    
##     0.941  a    
##     0.941  a    
##     0.899  a    
##     0.886  a    
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 12 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

# mc %>% write_sheet(ss = gs, sheet = "ModelClass")

plot2 <- mc %>% 
  plot_smr(type = "bar"
           , x = "architecture"
           , xlab = "Architecture" 
           , y = "emmean"
           , ylimits = c(0, 1.1, 0.2)
           , ylab = "Accuracy prediction"
           , group = "dense layers"
           , glab = "Dense layers"
           , sig = ".group"
           # , error = "SE"
           ) +
  facet_wrap(. ~ `activation function`, ncol = 3) +
  theme(legend.position = c(0.15, 0.97)
        , legend.direction = "horizontal"
        , strip.text = element_text(size = 12)
        , axis.text.x = element_text(size = 9)
        )

tab <-  mc %>% 
  merge(., mdl_class_info, by = c("architecture", "dense layers", "activation function"), all.x = )

# tab %>% 
#   mutate(across(where(is.numeric), ~round(., 3))) %>% 
#   select(model, everything()) %>% 
#   select(!c(df, lower.CL, upper.CL, num)) %>% 
#   rename(Model = "model"
#          , ste = "SE", sig = ".group", accuracy = "emmean") %>% 
#   mutate(across(sig, ~trimws(.))) %>% 
#   arrange(Model) %>% 
#   write_sheet(ss = gs, sheet = "tab3") 

plot <- plot_grid(plot1, plot2
                  , ncol = 1
                  , labels = "auto"
                  ) 

plot %>% 
  ggsave2(filename = "files/Fig2.jpg"
          , plot = .
          , dpi = 300
          , width = 10
          , height = 8
          )

include_graphics("files/Fig2.jpg")

9 Figure 2 (pipeline)

up <- 1

figure <- label_layout(size = c(20, 26), border_color = NULL) %>% 
  # include_text(value = "PhenomQuinoa", size = 20, position = c(10, 25.5)) %>% 
  include_text(value = "(a) Panicle images", size = 16, position = c(5, 24.5 + up)) %>% 
  include_image(value = "annotation/pipeline/orig_ATP.jpg"
                , position = c(3.5, 22 + up) , size = c(2.5, 4)
                ) %>% 
  include_image(value = "annotation/pipeline/orig_amaranth.jpg"
                , position = c(6.2, 22 + up) , size = c(3, 4)
                ) %>% 
  include_image(value = "annotation/pipeline/arrow.png"
                , position = c(10, 22 + up), size = c(1.2, 1.2)
                , opts = "image_rotate(90)*image_flop()"
                ) %>% 
  include_text(value = "(b) Image annotation", size = 16, position = c(15, 24.5 + up)) %>% 
  include_image(value = "annotation/pipeline/annotation.png"
                , position = c(15, 22 + up), size = c(5, 5)
                ) %>% 
  include_image(value = "annotation/pipeline/arrow.png"
                , position = c(15, 19 + up), size = c(1.2, 1.2)
                , opts = "image_rotate(0)*image_flop()"
                ) %>% 
  include_text(value = "(c) Segmentation model training:\n16 models", size = 16, position = c(15, 17.7 + up)) %>% 
  include_image(value = "annotation/pipeline/model.png"
                , position = c(15, 15 + up), size = c(6.5, 6.5)
                ) %>% 
  include_image(value = "annotation/pipeline/arrow.png"
                , position = c(10, 15 + up), size = c(1.2, 1.2)
                , opts = "image_rotate(-90)*image_flop()"
                ) %>% 
  include_text(value = "(d) Best segmentation model\n(segmentation-09)", size = 15, position = c(5, 17.7 + up)) %>% 
  include_image(value = "annotation/pipeline/mseg_ATP.jpg"
                , position = c(3.5, 15 + up)
                , size = c(2.5, 4)
                ) %>% 
  include_text(value = "Panicle traits:\n- Lenght\n- Width\n- Area\n- RGB values\n- Indices"
             , size = 12, position = c(5, 15 + up), opts = list(hjust = 0, vjust = 0)) %>% 
  include_image(value = "annotation/pipeline/arrow.png"
              , position = c(5, 12 + up)
              , size = c(1.2, 1.2), opts = "image_rotate(0)*image_flop()"
              ) %>% 
  include_text(value = "(e) Classification model training:\n12 models", size = 16, position = c(5, 10.7 + up)) %>% 
  include_image(value = "annotation/pipeline/model.png"
                , position = c(5, 7.7 + up), size = c(6.5, 6.5)
                ) %>% 
  include_image(value = "annotation/pipeline/arrow.png"
                , position = c(10, 7.7 + up), size = c(1.2, 1.2)
                , opts = "image_rotate(90)*image_flop()"
                ) %>% 
  include_text(value = "(f) Best classification model\n(classification-08)", size = 16, position = c(15, 10.8 + up)) %>% 
  include_image(value = "annotation/pipeline/mclass_ATP.jpg"
                , position = c(14, 7.8 + up)
                , size = c(2.5, 4)
                ) %>% 
  include_text(value = "Panicle shape:\n- Glomerulate\n- Amarantiform"
               , size = 12, position = c(15.5, 7.8 + up), opts = list(hjust = 0, vjust = 0)) %>% 
  include_text(value = "(g) Phenotyping\n(merge data)", size = 16, position = c(10, 6.5)) %>% 
  include_image(value = "annotation/pipeline/marrow.png"
            , position = c(10, 4.0)
            , size = c(3, 3), opts = 'image_transparent("white")'
            ) %>% 
  include_text(value = "(h) Analysis", size = 16, position = c(10, 2.3)) %>% 
  include_text(value = "- Genetic parameters\n- Heritability\n- QTL mapping*\n- GWAS*"
             , size = 12, position = c(8.5, 0.9), opts = list(hjust = 0, vjust = 0))
  
figure %>% label_print()

  
fig <- figure %>% 
  label_print(filename = "files/Fig2", mode = "c")

10 ModelSeg vs. ImageJ

imgj <- gs %>% 
  range_read("imagej") %>% 
  dplyr::mutate(across(everything(), as.character)) %>% 
  dplyr::mutate(rgb = dplyr::case_when(
    rgb == "1" ~ "R"
    , rgb == "2" ~ "G"
    , rgb  == "3" ~ "B"
  )) %>% 
  pivot_longer(!c(1:3)) %>% 
  mutate(name = case_when(
    name %in% "mean" ~ paste("panicle", rgb, name, sep = "_")
    , name %in% "stdev" ~ paste("panicle", rgb, name, sep = "_")
    , name %in% "mode" ~ paste("panicle", rgb, name, sep = "_")
    , TRUE ~ name
  )) %>% 
  mutate(value = case_when(
    rgb %in% "R" & name %in% c("panicle_area", "perimeter", "BX", "BY", "panicle_width", "panicle_length", "panicle_shape") ~ value
    , str_detect(name, "mean|stdev|mode") ~ value
  )) %>% 
  drop_na(value) %>% 
  dplyr::select(!rgb) %>% 
  pivot_wider() %>% 
  mutate(panicle_shape = case_when(
    panicle_shape %in% "amaranth" ~ 1
    , panicle_shape %in% "glomerulate" ~ 0
  )) %>%
  rename_with(~ gsub("panicle_", "imagej\n", .)) %>%
  dplyr::select(img_folder, img_name, starts_with("imagej"))
  
seg <- rs %>% 
  dplyr::select(img_folder, img_name, starts_with("panicle")) %>% 
  mutate(panicle_shape = case_when(
    panicle_shape %in% "amaranth" ~ 1
    , panicle_shape %in% "glomerulate" ~ 0
  )) %>% 
  rename_with(~ gsub("panicle_", "model\n", .))

comp <- merge(imgj, seg, all.x = T) %>% 
  mutate(across(everything(), as.character)) %>% 
  pivot_longer(!c(1:2)) %>% 
  separate(name, c("type", "trait"), sep = "\n") %>% 
  arrange(trait) %>% 
  unite(name, c(type, trait), sep = "\n") %>% 
  mutate(value = case_when(
    value %in% "glomerulate" ~ "0"
    , value %in% "amaranth" ~ "1"
    , TRUE ~ as.character(value)
  )) %>% 
  mutate(across( value, as.numeric)) %>% 
  pivot_wider() %>% 
  dplyr::select(!matches("number|_mode|detected|width.length|length.width"))

11 Trait correlation

library(ROCR)

predt <- comp %>% 
  rename_with(~ gsub("\n", " ", .)) %>% 
  select(contains("shape")) %>% 
  drop_na(`model shape`)

pred <- prediction(predictions = predt$`model shape`
                  , labels = predt$`imagej shape`)

perf1 <- performance(pred, "tpr", "fpr")

p1 <- ~ {
  plot(perf1,
     avg= "threshold",
     colorize=TRUE,
     lwd= 3)
}

perf2 <- performance(pred, "sens", "spec")

p2 <- ~{
  plot(perf2,
     avg= "threshold",
     colorize=TRUE,
     lwd= 3)
}

row2 <- list(p1, p2) %>% 
  plot_grid(plotlist = ., nrow = 1, labels = c("e", "f"))

library(caret)

pred <- table(predt$`model shape`, predt$`imagej shape`)

confusionMatrix(pred)
## Confusion Matrix and Statistics
## 
##    
##      0  1
##   0 29  0
##   1  1 12
##                                           
##                Accuracy : 0.9762          
##                  95% CI : (0.8743, 0.9994)
##     No Information Rate : 0.7143          
##     P-Value [Acc > NIR] : 0.00001297      
##                                           
##                   Kappa : 0.9431          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.9667          
##             Specificity : 1.0000          
##          Pos Pred Value : 1.0000          
##          Neg Pred Value : 0.9231          
##              Prevalence : 0.7143          
##          Detection Rate : 0.6905          
##    Detection Prevalence : 0.6905          
##       Balanced Accuracy : 0.9833          
##                                           
##        'Positive' Class : 0               
## 

library(ggside)

gside <- comp %>% 
  rename_with(~ gsub("\n", " ", .))

p1 <- gside %>% 
  inti::plot_raw(type = "s"
                 , x = "imagej length"
                 , y = "model length" 
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue") +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Image J - length (pixel)"
       , y = "Model - length (pixel)")

p2 <- gside %>% 
  inti::plot_raw(type = "s"
                 , x = "imagej width"
                 , y = "model width" 
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue", ) +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Image J - Width (pixel)"
       , y = "Model - Width (pixel)")

p3 <- gside %>% 
  inti::plot_raw(type = "s"
                 , x = "imagej area"
                 , y = "model area" 
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue", ) +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Image J - Area (pixel^2)"
       , y = "Model - Area (pixel^2)")


row1 <- list(p1, p2, p3) %>% 
  plot_grid(plotlist = ., nrow = 1, labels = "auto"
            , rel_widths = c(1, 1, 1.2))

list(row1, row2) %>% 
  plot_grid(plotlist = ., nrow = 2, rel_heights = c(1.2, 1)) %>% 
  ggsave2(filename = "files/Fig3.jpg"
          , width = 40
          , height = 25
          , units = "cm"
          )

list(row1, row2) %>% 
  plot_grid(plotlist = ., nrow = 2, rel_heights = c(1.2, 1)) %>% 
  ggsave2(filename = "files/Fig3.eps"
          , width = 40
          , height = 25
          , units = "cm"
          )

12 Multi-location trials analysis

rs$img_folder %>% unique()
## [1] "F6-camacani"  "F7-illpa"     "F8-camacani"  "F8-illpa"     "scale-images"

ilp2017 <- "https://docs.google.com/spreadsheets/d/1w7m54SKGGLEj1d4RB7JTDyS4957vZgzEg6XZREZRXGk/edit#gid=2102248480" %>% 
  as_sheets_id() %>% 
  range_read(sheet = "SC_org")

ilp17 <- ilp2017 %>% 
  dplyr::select(pht_flw, cross, line, r, sample,block, pnl_flw, pnd_flw
         , rdt_10
         ) %>% 
  unite("line", c(cross, line), sep = " ") %>% 
  mutate(img_name = paste0(pht_flw, ".jpg")) %>% 
  dplyr::select(!pht_flw)  %>% 
  mutate(nline = case_when(
    !str_detect(line, "x") ~ gsub(" .*", "\\1", line)
    , TRUE ~ line
  )) %>% 
  dplyr::select(!line) %>% 
  merge(., rs, sort = F, all.x = T) %>% 
  mutate(line = case_when(
    !is.na(nline) ~ nline
    , str_detect(img_name, pattern = "_") ~ sub("\\_.*", "", img_name)
  ), .after = img_name) 

#> illpa-2017-f7

ilp2017 <- ilp17 %>% 
  filter(img_folder %in% "F7-illpa") %>% 
  dplyr::select(line, r, block, sample
         , panicle_length.width, panicle_width.length, panicle_shape
         ) %>%
  mutate(across(c(line, r, block, sample), as.factor)) %>% 
  mutate(panicle_shape = case_when(
    panicle_shape %in% "glomerulate" ~ 0
    , panicle_shape %in% "amaranth" ~ 1
  )) %>% 
  mutate(across(where(is.character), as.numeric))  

model <- 5:length(ilp2017) %>% map(function(x) { 
  
  trait <- names(ilp2017)[x]
  
    H2cal(data = ilp2017
     , trait = trait
     , gen.name = "line"
     , fixed.model = "0 + (1|r) + (1|r:block) +  line"
     , random.model = "1 + (1|r) + (1|r:block) + (1|line)"
     , rep.n = 2
     , emmeans = F
     , plot_diag = T
     , outliers.rm = T
     , summary = T
     , trial = "F7-illpa"
     )
  
  })
## Linear mixed model fit by REML ['lmerMod']
## Formula: panicle_length.width ~ 1 + (1 | r) + (1 | r:block) + (1 | line)
##    Data: dt.rm
## Weights: weights
## 
## REML criterion at convergence: 5513.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.92990 -0.46067 -0.01391  0.50045  2.96695 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  line     (Intercept) 0.60427  0.7773  
##  r:block  (Intercept) 0.03282  0.1812  
##  r        (Intercept) 0.05658  0.2379  
##  Residual             0.26172  0.5116  
## Number of obs: 2760, groups:  line, 547; r:block, 20; r, 2
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   3.4917     0.1766   19.77

## Linear mixed model fit by REML ['lmerMod']
## Formula: panicle_width.length ~ 1 + (1 | r) + (1 | r:block) + (1 | line)
##    Data: dt.rm
## Weights: weights
## 
## REML criterion at convergence: -7794.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7438 -0.4408 -0.0306  0.4323  3.5783 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  line     (Intercept) 0.0056680 0.07529 
##  r:block  (Intercept) 0.0002224 0.01491 
##  r        (Intercept) 0.0003262 0.01806 
##  Residual             0.0019005 0.04359 
## Number of obs: 2712, groups:  line, 545; r:block, 20; r, 2
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.30583    0.01363   22.44

## Linear mixed model fit by REML ['lmerMod']
## Formula: panicle_shape ~ 1 + (1 | r) + (1 | r:block) + (1 | line)
##    Data: dt.rm
## Weights: weights
## 
## REML criterion at convergence: 2656.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5362 -0.4316 -0.1076  0.3806  2.5739 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  line     (Intercept) 0.154658 0.39327 
##  r:block  (Intercept) 0.006372 0.07982 
##  r        (Intercept) 0.015436 0.12424 
##  Residual             0.099681 0.31572 
## Number of obs: 2772, groups:  line, 548; r:block, 20; r, 2
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.44105    0.09151    4.82

tabsmr_ilp2017 <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
    }) %>% 
  bind_rows() %>% 
  select(!c(trial, rep, V.p, H2.s, H2.p, H2.c)) %>% 
  rename(Trait = trait) %>% 
  mutate(across(where(is.numeric), ~ round(., 3)))

# tabsmr_ilp2017 %>% sheet_write(ss = gs, sheet = "tabS1")

blues_ilp2017 <- 1:length(model) %>% map(function(x) {

      model[[x]][["blues"]] %>% dplyr::select(!matches("smith.w"))
    
      }) %>% 
  Reduce(function(...) merge(..., all = TRUE, by = c("trial", 'line')), .) %>% 
  bind_rows()

#> F8-camacani

cmc2018 <- rs %>% 
  filter(img_folder %in% "F8-camacani") %>% 
  separate(img_name, into = c("line", "suffix"), sep = "_", remove = FALSE, extra = "drop", fill = "right") %>% 
  dplyr::select(line
                , panicle_length.width, panicle_width.length
                , panicle_shape
                ) %>%
  mutate(across(c(line), as.factor)) %>% 
  mutate(panicle_shape = case_when(
    panicle_shape %in% "glomerulate" ~ 0
    , panicle_shape %in% "amaranth" ~ 1
  )) %>% 
  mutate(across(where(is.character), as.numeric)) %>% 
  mutate(trial = "F8-camacani", .before = line) 
  

#> F8-illpa

ilp2018 <- rs %>% 
  filter(img_folder %in% "F8-illpa") %>% 
  separate(img_name, into = c("line", "suffix"), sep = "_", remove = FALSE, extra = "drop", fill = "right") %>% 
  dplyr::select(line
                , panicle_length.width, panicle_width.length
                , panicle_shape
                ) %>%
  mutate(across(c(line), as.factor)) %>% 
  mutate(panicle_shape = case_when(
    panicle_shape %in% "glomerulate" ~ 0
    , panicle_shape %in% "amaranth" ~ 1
  )) %>% 
  mutate(across(where(is.character), as.numeric)) %>% 
  mutate(trial = "F8-illpa", .before = line) 

lines <- ilp2018 %>% select(line) %>% unique() %>% deframe() %>% as.vector()

# MLT data

mlt <- list(blues_ilp2017
            , cmc2018
            , ilp2018
            ) %>% 
  bind_rows() %>% 
  separate(trial, c("year", "location"), remove = F) %>%  
  mutate(across(where(is.character), as.factor)) 

mlt %>% names()
## [1] "trial"                "year"                 "location"            
## [4] "line"                 "panicle_length.width" "panicle_width.length"
## [7] "panicle_shape"

12.1 One stage analysis

mltbn <- list(ilp2017 %>% mutate(trial = "F7-illpa")
              , ilp2018
              , cmc2018 
              ) %>% 
  bind_rows() %>% 
  select(trial, everything()) 

model <- 6:length(mltbn) %>% map(function(x) { 
  
  trait <- names(mltbn)[x]
  
  mltbn %>% 
    H2cal(data = .
     , trait = trait
     , gen.name = "line"
     , env.name = "trial"
     , env.n = 3
     , fixed.model = "0 + (1|trial) + (1|line:trial) + line"
     , random.model = "1 + (1|trial) + (1|line:trial) + (1|line)"
     , rep.n = 2
     , emmeans = F
     , plot_diag = T
     , outliers.rm = switch(trait, panicle_shape = FALSE, TRUE)
     , summary = F
     )
  
  })

#> model parameter

model_tabsmr <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
    }) %>% bind_rows() %>% 
  select(!c(rep, V.p, repeatability, H2.p)) %>% 
  rename(Trait = trait) %>% 
  mutate(across(where(is.numeric), ~ round(., 3)))

model_tabsmr %>% web_table()

# model_tabsmr %>% sheet_write(ss = gs, sheet = "tabS2")

model_blups <- 1:length(model) %>% map(function(x) {

      model[[x]][["blups"]] %>% select(!matches("smith.w"))
    
      }) %>% Reduce(function(...) merge(..., all = TRUE, by = c('line')), .) %>% 
  bind_rows() %>% 
  pivot_longer(!line) %>% 
  pivot_wider()

model_blups %>% web_table(caption = "BLUPs")

model_blues <- 1:length(model) %>% map(function(x) {

      model[[x]][["blues"]] %>% select(1:2)
    
      }) %>% Reduce(function(...) merge(..., all = TRUE, by = c('line')), .) %>% 
  bind_rows() 

model_blues %>% web_table(caption = "BLUEs")

#>


plotvc <- model_tabsmr %>% 
  select(Trait, matches("V.g|V.e|H2")) %>% 
  select(where(~ any(. != 0))) %>% 
  pivot_longer(!c(Trait, H2.c)) %>% 
  group_by(Trait) %>% 
  mutate(percent = value/sum(value)*100) %>% 
  ungroup() %>% 
  mutate(H2 = round(H2.c, 2)) %>% 
  mutate(sig = paste0("H2 = ", round(H2.c, 2))) %>% 
  mutate(Trait = trimws(Trait)) %>% 
  mutate(Trait = gsub("\\.", "/", Trait)) %>% 
  mutate(Trait = gsub("_", "\n", Trait)) 

f4c <-  plotvc %>% 
    plot_smr(x = "Trait"
           , y = "percent"
           , group = "name"
           , xlab = ""
           , ylab = "Percentage ('%')"
           , glab = "Variance component"
           , ylimits = c(0, 109, 10)
           # , xrotation = c(45, 1, 1)
           ) +
  geom_bar(position="stack", stat="identity") +
  geom_text(data = plotvc, aes(x = Trait, label = sig, y = 105, fill=NULL)) 


lbls <- as_labeller(c("panicle_length.width" = "Panicle length/width"
                    , "panicle_width.length" = "Panicle width/length"
                    , "panicle_shape" = "Panicle shape"
                    ))

f4d <- model_blues %>%
  pivot_longer(!c(line)) %>% 
  ggplot(aes(x = value)) + 
    geom_histogram(aes(y=..density..), 
                   colour="black", fill="white") +
  scale_y_continuous(expand = c(0, 0)) + 
  geom_density(alpha=.1, fill="#FF6666") +
  facet_wrap(vars(name), scales = "free", nrow = 1
             , strip.position = "bottom"
             , labeller = lbls
             ) +
  theme_bw() +
  labs(x = "", y = "") 

13 Binomial distribution

library(lme4)
library(lmerTest)

bmodel1 <- lmer(panicle_shape ~ 1 + (1|trial) + (1|line:trial) + (1|line)
                , data = mltbn)

bmodel2 <- glmer(panicle_shape ~ 1 + (1|trial) + (1|line:trial) + (1|line)
                , family = binomial(link = "logit")
                , data = mltbn)

AIC(bmodel1, bmodel2)
##         df      AIC
## bmodel1  5 2909.387
## bmodel2  4 2807.708

#> lmer

g.ran <- bmodel1
gen.name <- "line"
env.name <- "trial"
env <- 3
rep <- 2

g.ran %>% summary()
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: panicle_shape ~ 1 + (1 | trial) + (1 | line:trial) + (1 | line)
##    Data: mltbn
## 
## REML criterion at convergence: 2899.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6244 -0.2477 -0.1381  0.3732  2.6416 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  line:trial (Intercept) 0.06386  0.2527  
##  line       (Intercept) 0.07809  0.2794  
##  trial      (Intercept) 0.01098  0.1048  
##  Residual               0.11118  0.3334  
## Number of obs: 2821, groups:  line:trial, 597; line, 549; trial, 3
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)  0.33082    0.07479 1.72403   4.424   0.0618 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vc.g <- c(lme4::VarCorr(g.ran)[[gen.name]])

gxl <- paste(gen.name, env.name, sep = ":")

vc.gxl <- g.ran %>%
      lme4::VarCorr() %>%
      tibble::as_tibble() %>%
      dplyr::filter(grp == gxl) %>%
      dplyr::pull(vcov)

vc.e <- g.ran %>%
    lme4::VarCorr() %>%
    tibble::as_tibble() %>%
    dplyr::filter(grp == "Residual") %>%
    dplyr::pull(vcov)

vc.p = (vc.g + vc.gxl/env + vc.e/(env*rep))

H2.s = vc.g/vc.p

H2.s 
## [1] 0.6622923

#> lmer = 0.66

#> glmer

g.ran <- bmodel2
gen.name <- "line"
env.name <- "trial"
env <- 3
rep <- 2

g.ran %>% summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: panicle_shape ~ 1 + (1 | trial) + (1 | line:trial) + (1 | line)
##    Data: mltbn
## 
##      AIC      BIC   logLik deviance df.resid 
##   2807.7   2831.5  -1399.9   2799.7     2817 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7231 -0.1930 -0.1523  0.2921  2.8102 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  line:trial (Intercept) 9.61     3.10    
##  line       (Intercept) 7.13     2.67    
##  trial      (Intercept) 0.00     0.00    
## Number of obs: 2821, groups:  line:trial, 597; line, 549; trial, 3
## 
## Fixed effects:
##             Estimate Std. Error z value      Pr(>|z|)    
## (Intercept)  -1.4855     0.2528  -5.876 0.00000000419 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

vc.g <- c(lme4::VarCorr(g.ran)[[gen.name]])

gxl <- paste(gen.name, env.name, sep = ":")

vc.gxl <- g.ran %>%
      lme4::VarCorr() %>%
      tibble::as_tibble() %>%
      dplyr::filter(grp == gxl) %>%
      dplyr::pull(vcov)

vc.e <- g.ran %>%
    lme4::VarCorr() %>%
    tibble::as_tibble() %>%
    dplyr::filter(grp == "Residual") %>%
    dplyr::pull(vcov)

vc.e <- 0

vc.p = (vc.g + vc.gxl/env + vc.e/(env*rep))

H2.s = vc.g/vc.p

H2.s 
## [1] 0.6899932

#> glmer = 0.69

13.1 Panicle distribution by type

pdt <- list(ilp2017 %>% mutate(trial = "F7-illpa", .before = line)
            , cmc2018
            , ilp2018
            ) %>% 
  bind_rows() %>% 
  mutate(panicle_shape = case_when(
    panicle_shape %in% 0 ~ "glomerulate"
    , panicle_shape %in% 1 ~ "amaranth"
  )) %>% 
  separate(trial, c("generation", "location"), remove = F) 
  
panper <- pdt %>% 
  group_by(generation, location, panicle_shape) %>% 
  summarise(n = n()) %>% 
  mutate(per = (n/ sum(n)) * 100) %>% 
  mutate(across(per, ~round(., 1)))

# panper %>% sheet_write(ss = gs, sheet = "tab4")

#> plot

pie <- panper %>% 
  mutate(csum = rev(cumsum(rev(per))), 
         pos = per/2 + lead(csum, 1),
         pos = if_else(is.na(pos), per/2, pos))

f4a <- pie %>% 
  ggplot(aes(x = "" , y = per, fill = panicle_shape)) +
  geom_col(width = 1, color = 1) +
  coord_polar(theta = "y") +
  scale_fill_brewer(palette = "Pastel2") +
  ggrepel::geom_label_repel(data = pie, 
                   aes(y = pos, label = paste0(per, "% (n =", n , ")")),
                   size = 2.5, nudge_x = 1, show.legend = FALSE) +
  guides(fill = guide_legend(title = "Panicle Shape", title.position = "top")) +
  facet_wrap(generation ~ location, ncol = 1) +
  theme_void() +
  theme(legend.position="bottom") 

f4a

#> 

modelpan <- pdt %>% 
  aov(panicle_length.width ~ trial*panicle_shape, data = .)

Anova(modelpan,type="III", singular.ok = TRUE)
## Anova Table (Type III tests)
## 
## Response: panicle_length.width
##                      Sum Sq   Df   F value                Pr(>F)    
## (Intercept)          834.12    1 1347.7715 < 0.00000000000000022 ***
## trial                  2.69    2    2.1740                0.1139    
## panicle_shape         15.65    1   25.2905          0.0000005241 ***
## trial:panicle_shape    0.10    2    0.0838                0.9196    
## Residuals           1742.18 2815                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res1 <- emmeans::emmeans(modelpan, ~ panicle_shape) %>% 
  multcomp::cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, ~trimws(.)))
  
f4b <- pdt %>% 
  plot_raw(x = "panicle_shape"
           , y = "panicle_length.width"
           , legend = "none"
           , ylab = "Panicle length/width"
           , xlab = "Panicle shape"
           ) +
  annotate("text", y = 7.5
           , x = res1$panicle_shape
           , label = res1$.group)

13.2 model selection

#> model comparison

m1 <- lme4::lmer(panicle_length.width ~ 1 + (1|location) + (1|year) + (1|line:location) + (1|line:year) + (1|line:year:location) + (1|line), data = mlt)

m2 <- lme4::lmer(panicle_length.width ~ 1 + (1|trial) + (1|line:trial) + (1|line), data = mlt)

aic.1 <- tibble(model = c("m1", "m2"), model_formula = c(
  "trait ~ 1 + (1|location) + (1|year) + (1|line:location) + (1|line:year) + (1|line:year:location) + (1|line)"
  , "trait ~ 1 + (1|trial) + (1|line:trial) + (1|line)"
)) 

#> error for low level in groups

m3 <- lme4::lmer(panicle_length.width ~ 1 + (1|location) + (1|year) + (1|line:location) + (1|line:year)  + (1|line), data = mlt)

m4 <- lme4::lmer(panicle_length.width ~ 1 + (1|trial) + (1|line), data = mlt)

aic.2 <- AIC(m3, m4) %>% 
  data.frame() %>% 
  rownames_to_column("model") %>% 
  mutate(model_formula = case_when(
    model == "m3" ~ "trait ~ 1 + (1|location) + (1|year) + (1|line:location) + (1|line:year)  + (1|line)"
    , model == "m4" ~ "trait ~ 1 + (1|trial) + (1|line)"
  ))

model_aic <- bind_rows(aic.1, aic.2) %>% 
  select(model, df, AIC, everything())

# model_aic %>% sheet_write(ss = gs, sheet = "STab1")

model_aic %>% web_table()

13.3 quantitative-gentics parameters

#> Complete model

model <- 5:length(mlt) %>% map(function(x) { 
  
  trait <- names(mlt)[x]
  
  mlt %>% 
    H2cal(data = .
     , trait = trait
     , gen.name = "line"
     , env.name = "location"
     , env.n = 2
     , year.name = "year"
     , year.n = 2
     , fixed.model = "0 + (1|location) + (1|year) + (1|line:location) + (1|line:year) + line"
     , random.model = "1 + (1|location) + (1|year) + (1|line:location) + (1|line:year) + (1|line)"
     , rep.n = 2
     , emmeans = F
     , plot_diag = T
     , outliers.rm = F
     , summary = F
     )
  
  })

#> model parameter

model_tabsmr <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
    }) %>% bind_rows() %>% 
  select(!c(rep, V.p, repeatability, H2.p, V.gxy, V.gxlxy)) %>% 
  rename(Trait = trait) %>% 
  mutate(across(where(is.numeric), ~ round(., 3)))

model_tabsmr %>% web_table()

# model_tabsmr %>% sheet_write(ss = gs, sheet = "STable7")

model_blups <- 1:length(model) %>% map(function(x) {

      model[[x]][["blups"]] %>% select(!matches("smith.w"))
    
      }) %>% Reduce(function(...) merge(..., all = TRUE, by = c('line')), .) %>% 
  bind_rows() %>% 
  pivot_longer(!line) %>% 
  pivot_wider()

model_blups %>% web_table(caption = "BLUPs")

model_blues <- 1:length(model) %>% map(function(x) {

      model[[x]][["blues"]] %>% select(1:2)
    
      }) %>% Reduce(function(...) merge(..., all = TRUE, by = c('line')), .) %>% 
  bind_rows() 

model_blues %>% web_table(caption = "BLUEs")

13.4 variance component plot

plotvc <- model_tabsmr %>% 
  select(Trait, matches("V.g|V.e|H2")) %>% 
  select(where(~ any(. != 0))) %>% 
  pivot_longer(!c(Trait, H2.c)) %>% 
  group_by(Trait) %>% 
  mutate(percent = value/sum(value)*100) %>% 
  ungroup() %>% 
  mutate(H2 = round(H2.c, 2)) %>% 
  mutate(sig = paste0("H2 = ", round(H2.c, 2))) %>% 
  mutate(Trait = trimws(Trait)) %>% 
  mutate(Trait = gsub("\\.", "/", Trait)) %>% 
  mutate(Trait = gsub("_", "\n", Trait)) 

f4c <-  plotvc %>% 
    plot_smr(x = "Trait"
           , y = "percent"
           , group = "name"
           , xlab = ""
           , ylab = "Percentage ('%')"
           , glab = "Variance component"
           , ylimits = c(0, 109, 10)
           # , xrotation = c(45, 1, 1)
           ) +
  geom_bar(position="stack", stat="identity") +
  geom_text(data = plotvc, aes(x = Trait, label = sig, y = 105, fill=NULL)) 

13.5 Trait distribution

lbls <- as_labeller(c("panicle_length.width" = "Panicle length/width"
                    , "panicle_width.length" = "Panicle width/length"
                    , "panicle_shape" = "Panicle shape"
                    ))

f4d <- model_blups %>%
  pivot_longer(!c(line)) %>% 
  ggplot(aes(x = value)) + 
    geom_histogram(aes(y=..density..), 
                   colour="black", fill="white") +
  scale_y_continuous(expand = c(0, 0)) + 
  geom_density(alpha=.1, fill="#FF6666") +
  facet_wrap(vars(name), scales = "free", nrow = 1
             , strip.position = "bottom"
             , labeller = lbls
             ) +
  theme_bw() +
  labs(x = "", y = "") 

14 Figure 4

f4bc <- plot_grid(f4b
                 , f4c
                 , labels = c("b", "c")) 


f4bcd <- plot_grid(f4bc, f4d
                , nrow = 2
                , labels = c("", "d")
                )

f4 <- plot_grid(f4a, f4bcd
                , nrow = 1
                , labels = c("a", "")
                , rel_widths = c(1, 4)
                )

f4 %>% 
  ggsave(plot = ., "files/Fig4.jpg", width = 25, height = 20, units = "cm")

f4 %>% 
  ggsave(plot = ., "files/Fig4.eps", width = 25, height = 20, units = "cm")

15 Figure 5

segscale <- seg %>% 
  filter(img_folder == "scale-images")

scaledt <- gs %>% 
  range_read("scale") %>% 
  filter(!id %in% 236) %>% 
  pivot_longer(!c(1:7)) %>% 
  filter(grepl("white", name)) %>%
  drop_na("panicle-shape") %>% 
  mutate(`panicle-shape` = case_when(
    `panicle-shape` %in% c("2", "3") ~ "1"
    , `panicle-shape` == c("1") ~ "0"
  )) %>% 
  mutate(across(`panicle-shape`, as.numeric)) %>% 
  rename(img_name = "value") %>% 
  mutate(`panicle-width.length` = `panicle-width`/`panicle-length`) %>% 
  mutate(`panicle-length.width` = `panicle-length`/ `panicle-width`) %>% 
  rename_with(~gsub("panicle-", "scale\n", .))

dtscale <- merge(scaledt, segscale, by.x = "img_name", all.x = T) %>% 
  dplyr::select(matches("img_name|length|width|shape")) %>% 
  mutate(across(everything(), as.character)) %>% 
  drop_na(`model\nshape`) %>% 
  pivot_longer(!img_name) %>% 
  separate(name, c("model", "trait"), sep = "\n") %>% 
  arrange(trait) %>% 
  unite("name", c("model", "trait")) %>% 
  pivot_wider() %>% 
  mutate(across(!img_name, ~as.numeric(.))) 

library(ggside)

p5b <- dtscale %>% 
  inti::plot_raw(type = "s"
                 , x = "scale_length" 
                 , y = "model_length"
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue") +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Real - length (cm)"
       , y = "Model - length (pixel)")

p5c <- dtscale %>% 
  inti::plot_raw(type = "s"
                 , x = "scale_width" 
                 , y = "model_width"
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue") +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Real - width (cm)"
       , y = "Model - width (pixel)")

p5d <- dtscale %>% 
  inti::plot_raw(type = "s"
                 , x = "scale_width.length" 
                 , y = "model_width.length"
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue") +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Real - widwidth/lengthth (cm)"
       , y = "Model - width/length (pixel)")

p5e <- dtscale %>% 
  inti::plot_raw(type = "s"
                 , x = "scale_length.width" 
                 , y = "model_length.width"
                 ) +
  geom_xsidedensity(alpha = 0.25, fill = "red", outline.type = "lower") +
  geom_ysidedensity(alpha = 0.25, fill = "blue") +
  scale_ysidex_continuous(guide = guide_axis(angle = 90), position = "top", n.breaks = 3) +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("eq", "R2"))) +
  theme(ggside.panel.scale.x = 0.25
        , ggside.panel.scale.y = 0.25
        ) +
  labs(x = "Real - length/width (cm)"
       , y = "Model - length/width (pixel)")


library(ROCR)

pred <- prediction(predictions = dtscale$model_shape
                  , labels = dtscale$scale_shape)

perf1 <- performance(pred, "tpr", "fpr")


pred <- table(dtscale$scale_shape, dtscale$model_shape)

confusionMatrix(pred)
## Confusion Matrix and Statistics
## 
##    
##      0  1
##   0 36  1
##   1  0 16
##                                           
##                Accuracy : 0.9811          
##                  95% CI : (0.8993, 0.9995)
##     No Information Rate : 0.6792          
##     P-Value [Acc > NIR] : 0.00000003257   
##                                           
##                   Kappa : 0.956           
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.9412          
##          Pos Pred Value : 0.9730          
##          Neg Pred Value : 1.0000          
##              Prevalence : 0.6792          
##          Detection Rate : 0.6792          
##    Detection Prevalence : 0.6981          
##       Balanced Accuracy : 0.9706          
##                                           
##        'Positive' Class : 0               
## 

p5f <- ~ {
  plot(perf1,
     avg= "threshold",
     colorize=TRUE,
     lwd= 3)
}

p5a <- "annotation/figures/fig5a.jpg" %>% 
  jpeg::readJPEG() %>% 
  grid::rasterGrob()

p5ad <- list(p5b, p5c, p5d, p5e) %>% 
  plot_grid(plotlist = ., ncol = 2, labels = "auto")

p5fe <- list(p5f, p5a) %>% 
  plot_grid(plotlist = ., labels = c("e", "f"), ncol = 1
              , rel_heights = c(0.7, 1)) 

list(p5ad, p5fe) %>% 
  plot_grid(plotlist = ., ncol = 2, rel_widths = c(2, 1)) %>% 
  ggsave(plot = .
         , "files/Fig5.jpg", width = 30, height = 20, units = "cm")

list(p5ad, p5fe) %>% 
  plot_grid(plotlist = ., ncol = 2, rel_widths = c(2, 1)) %>% 
  ggsave(plot = .
         , "files/Fig5.eps", width = 30, height = 20, units = "cm")

16 Supplementary Figure 2

dir <- "annotation/imagej/" %>% 
  list.files(pattern = ".png", full.names = T)
  
plots <- 1:length(dir) %>% map( function(x) {
  
  dir[[x]] %>% 
    image_import() %>% 
    grid::rasterGrob()
  
})

fig <- cowplot::plot_grid(plotlist = plots
                   , ncol = 2
                   , labels = c("(a) Parameters selection"
                                , "(b) Image analysis")
                   , hjust = 0
                   , label_y = 1.07
                   , label_x = 0.01
                   , label_size = 25
                   ) +
  theme(plot.margin = unit(c(1, 0, 0, 0), "cm"))

fig %>% 
  ggsave2(filename = "files/FigS2.jpg"
          , width = 20
          , height = 7.5
          )

"files/FigS2.jpg" %>% include_graphics()

17 Supplementary Figure 3

dir <- "annotation/figures/" %>% 
  list.files(pattern = "FigS2a|b", full.names = T)
  
plots <- 1:length(dir) %>% map( function(x) {
  
  dir[[x]] %>% 
    image_import() %>% 
    grid::rasterGrob()
  
})

fig <- cowplot::plot_grid(plotlist = plots
                   , ncol = 2
                   , labels = "auto"
                   , label_size = 20
                   ) 

fig %>% 
  ggsave2(filename = "files/FigS3.jpg"
          , width = 12
          , height = 8
          )

"files/FigS3.jpg" %>% include_graphics()

18 Image pdf to png

pdf2png <- list.files("files/", pattern = "Fig.*pdf", full.names = T) 

imgs <- 1:length(pdf2png) %>% map(\(x) {
  
  filename <- pdf2png[x] %>% gsub(".pdf", "\\1.png", .) 
  
  pdf2png[x] %>% 
    image_read_pdf() %>% 
    image_write(format = "png", filename)
  
})